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Abstract: In this paper we investigate some stochastic models for 
tumor-immune systems. To describe these models, we used a Wiener 
process, as the noise has a stabilization effect. Their dynamics are studied 
in terms of stochastic stability in the equilibrium points, by constructing 
the Lyapunov exponent, depending on the parameters that describe the 
model. Stochastic stability was also proved by constructing a Lyapunov 
function. We have studied and and analyzed a Kuznetsov- Taylor like 
stochastic model and a Bell stochastic model for tumor-immune systems. 
These stochastic models are studied from stability point of view and they 
were represented using the second Euler scheme and Maple 12 software. 



1 Introduction 

Stochastic modeling plays an important role in many branches of 
science. In many practical situations perturbations appear and these are 
expressed using white noise, modeled by Brownian motion. We will study 
stochastic dynamical systems that are used in medicine, in describing a 
tumor behavior, but still we don't know much about the mechanism 
of destruction and establishment of a cancer tumor, because a patient 
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may experience tumor regression and later a relapse can occur. The 
need to address not only preventative measures, but also more successful 
treatment strategies is clear. Efforts along these lines are now being 
investigated through immunotherapy ([6], [21], [23]). 

This tumor-immune study, from theoretical point of view, has been 
done for two cell populations: effector cells and tumor cells. It was 
predicted a threshold above which there is uncontrollable tumor growth, 
and below which the disease is attenuated with periodic exacerbations 
occurring every 3-4 months. There was also shown that the model does 
have stable spirals, but the Dulac-Bendixson criterion demonstrates that 
there are no stable closed orbits. It is consider ODE's for the populations 
of immune and tumor cells and it is shown that survival increases if the 
immune system is stimulated, but in some cases an increase in effector 
cells increases the chance of tumor survival. 

In the last years, stochastic growth models for cancer cells were de- 
veloped, we mention the papers of W.Y. Tan and CW. Chen [20], N. 
Komarova, G. Albano and V.Giorno [2], L. Ferrante, S. Bompadre, L. 
Possati and L. Leone [7J, A. Boondirek Y. Lenbury, J. Wong-Ekkabut, 
W. Triampo, I.M. Tang, P. Picha [1]. 

Our goal in this paper is to construct stochastic models and to analyze 
their behavior around the equilibrium point. In these points stability is 
studied by analyzing the Lyapunov exponent, depending of the parame- 
ters of the models. Numerical simulations are done using a deterministic 
algorithm with an ergodic invariant measure. In this paper the authors 
studied and analyzed two stochastic models. In Section 2, we considered 
a Kuznetsov and Taylor stochastic model. Beginning from the classical 
one, we have studied the case of positive immune response. We gave the 
stochastic model and we analyzed it in the equilibrium points. Numerical 
simulations for this new model are presented in Section 2.1. In Section 3 
we presented a general family of tumor-immune stochastic systems and 
from this general representation we analyzed Bell model. We wrote this 
model as a stochastic model, using Annexe 1, and we discussed its behav- 
ior around the equilibrium points. We have proved stochastic stability 
around equilibrium point using two methods. The first one consists of ex- 
pressing the Lyapunov exponent, and then drawing the conclusion when 
the considered system is stable. The second method is a way of construct- 
ing a Lyapunov function and determining sufficient conditions such that 
the system is stable. Numerical simulations were done using the soft- 
ware Maple 12 and we implemented the second order Euler scheme for a 
representation of the discussed stochastic models. 
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2 Kuznetsov and Taylor stochastic model 



The study of tumor- immune interaction is determined by the behav- 
ior of two populations of cells: effector cells and tumor cells. We will 
construct the stochastic models using well known deterministic models 
and we analyze stochastic stability around the equihbrium points. The 
analysis is done using Lyapunov exponent method. 

We will begin our study from the deterministic model of Kuznetsov 
and Taylor [Tl]. This model describes the response of effector cells to 
the growth of tumor cells and takes into consideration the penetration of 
tumor cells by effector cells, that causes the interaction of effector cells. 
This model can be represented in the following way: 

f x{t) = ai - a2x{t) + a3x{t)y{t), . . 

\ y{t) = biymi-hy{t))-x{t)y{t), ^"-^ 

where initial conditions are x(0) = Xq > 0, y{0) = yo > and 03 is the 
immune response to the appearance of the tumor cells. 

In this paper we consider the case of 03 > 0, that means that immune 
response is positive. For the equilibrium states Pi and P2, we study the 
asymptotic behavior with respect to the parameter oi in ([T]). For bia2 < 
Oi, the system ([T]) has the equilibrium states Pi{xi,yi) and ^2(2^25 1/2)5 
with 

xi = -, z/i = 0, (2) 

X2 = (6i(a3-&2a2) + y(A))/(2a3), 1/2 = (61(03 + &2a2) - y(A)/(26i62a3) 

(3) 

where A = 6^(6202 — ^3)^ + 46i620ia3- 

We associate a stochastic system of differential equations to the ordi- 
nary system of differential equations ([T]). 

In [13] it is shown that there is an aio, such that if < ai < aio, then 
the equilibrium state Pi is asymptotically stable, and for ai > aiQ the 
equilibrium state Pi is unstable. If ai < aio, then the equilibrium state 
P2 is unstable and for ai > oio it is asymptotically stable. 

Let us consider (f2, 3^t>o, 7) a filtered probability space and {W{t))t>o 
a standard Wiener process adapted to the filtration {3^)t>o- Let {X{t, u) = 
{x(t),y{t))}t>Q be a stochastic process. 

The system of Ito equations associated to system ([T]) is given by 

f x{t)=xo + Jl^{ai - a2x{s) + a3x{s)y{s))ds + J^^ gi{x{s),y{s))dW{s), 
I y{t) = yo + /q((6i?/(s)(1 - b2y{s)) - x{s)y{s))ds + J^^ g2{x{s),y{s))dW{s), 

(4) 
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where the first integral is a Riemann integral, and the second one is an 
Ito integral. {W(t)}t>o is a Wiener process [T7] . 

The functions gi{x{t),y{t)) and g2{x{t),y{t)) are given in the case 
when we are working in the equilibrium state. In Pi those functions have 
the following form 

gi{x(t), y{t)) = biixit) + buyit) + cn, 

(5) 

92{x{t), y{t)) = b2ix{t) + b22y{t) + C21, 

where 

cii = -biiXi - 6i2?/i, C21 = -b2ixi - 622^1- (6) 

In the equilibrium state P2, the functions gi{x{t),y{t)) and g2{x(t),y(t)) 
are given by 

gi{x{t), y{t)) = biix{t) + bi2y{t) + C12, 

(7) 

g2{,x(t), y{t)) = b2ix(t) + b22y{t) + C22, 

where 

C12 = -biiX2 - fei2l/2, C22 = -6213:2 - b22y2- (8) 

The functions gi{x{t),y{t)) and g2{x{t),y{t)) represent the volatilisa- 
tions of the stochastic equations and they are the therapy test functions. 

2.1 The analysis of SDE ([41). Numerical simulation. 

Using the formulae from Annexe 1, Annexe 2, and Maple 12 soft- 
ware, we get the following results, illustrated in the below figures. For 
numerical simulations we use the following values of parameters: 

ai = 0.1181, a2 = 0.3747, 03 = 0.01184, 61 = 1.636, 62 = 0.002. 

The matrices A and B are given, in the equilibrium point Pi(^, 0) by 

'-02 + 032/1 a^iXi \ ^ flO -2 



-yi bi - 2622/1 - ' V 2 10 

In a similar way, matrices A and B are defined in the other equilibrium 
point 

^ / (-fei(&2Q2 -Q3) + VA) (&i(&2a2 + Q3) - VA )- 
203 ' 2616203 

with A = 61(6202 — a^Y + 461620103. 

Using second order Euler scheme, for the ODE system ([T]) and SDE 
system we get the following orbits presented in the figures above. 



4 




200 400 600 800 1 ,000 1 00 200 300 400 500 



Figure 1: {n,x{n)) in Pi for ODE ^ Figure 2: {n,x{n,uj)) in Pi for SDE Q 
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Figure 3: {n,y{n)) in Pi for ODE ^ Figure 4: {n,y{n,uo)) in Pi for SDE (g]) 




Figure 5: {x{n),y{n)) in Pi for ODE ([T]) Figure 6: {x{n,uj),y{n,uj)) in Pi for SDE 
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Figure 7: {n,x{n)) in P2 for ODE ^ Figure 8: {n,x{n,uj)) in P2 for SDE (g]) 
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Figure 9: {n,y{n)) in P2 for ODE ^ Figure 10: {n,y{n,Lu)) in P2 for SDE (g]) 




Figure 11: {x{n),y{n)) in P2 for ODE ((T|) Figure 12: (a;(n, oj), y(n, tj)) in P2 for SDE 
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The Lyapunov exponent variation, with the variable parameter bn = a, 
is given in Figure 13 for Pi, and in Figure 14 for the equihbrium point 

P2. 




Figure 13: {a, A(a)) in Pi for ODE Figure 14: {a, \{a)) in P2 for SDE 



The Lyapunov exponent, for the equihbrium point Pi is negative, so 
Pi is asymptotically stable for each a G M. For the second equilibrium 
point P2, this point is asymptotically stable for all values when A < 0, 
that means that P2 is unstable for all a e {—00, —1.8) U (1.8, 00). 

3 A general family of tumor- immune stochas- 
tic systems 

A Volterra-like model was proposed in [18] for the interaction be- 
tween a population of tumor cells (whose number is denoted by x) and 
a population of lymphocyte cells {y) and it is given by 

r x{t) = ax{t) - hx{t)y{t), , . 

\ y{t) = dx{t)y{t) - fx{t) - kx{t) + u, ^ ^ 

where the tumor cells are supposed to be in exponential growth (which is, 
however, a good approximation only for the initial phases of the growth) 
and the presence of tumor cells implies a decrease of the "input rate" of 
lymphocytes. 

A general representation for such models can be considered in the 
form given by d'Onofrio in [B]: 

r x{t) = fi{x{t),y{t)), y{t) = f2{x{t),y{t)), 

\ x(0) = xo, 2/(0) = 2/0, ^ ' 
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where x is the number of tumor cells, y the number of effector cells of 
immune system and 

fi{x,y) = x{hix - h2xy), 

(11) 

/2(x, y) = {hsx - hix)y + h^x. 

The functions hi,h2,h3, /14, hr, are given such that the system (ITOl) admits 
the equilibrium point Pi{xi, yi), with xi = 0, yi > 0, and ^2(3^2, 2/2)5 with 
xi 7^ 0, 2/2 > 0. 

Deterministic models of this general form are the following 

Volterra model |22] if hi{x) = ai, h2{x) = a2X, h^{x) = b^x, h^^i^x) = 
62 and h^i^x) = —bix; 

Bell model hi{x) = aix, h2{x) = a2X, h^{x) = bix, /i4(x) = 63 and 
hix) = -b2X + 64; 

Stepanova model [19] with/;,i(a;) = ai, h2{x) = 1, h^lx) = bix, h^lx) = 
b and h^i^x) = —b2X + 64; 

Vladar-Gonzalez model [21] if in (ITOl) we consider hi{x) = \og{K/x), 
h2{x) = 1, hsi^x) = bix, h^{x) = 62 + b^x"^ and h^{x) = 1; 

Exponential model | |23] if in ( ITOl) we consider hi{x) = 1, h2{x) = 1, 
hsi^x) = bix, hi{x) = 62 + b'ix'^, and h^{x) = 1; 

Logistic model | |15] if in (ITOl) we consider hi{x) = 1 — — , h2{x) = 1, 
hsix) = bix, hi{x) = 62 + ^33;^, and h^{x) = 1. 

The analysis of these models was proven also using numerical simulations. 

For a considered filtered probability space (f2, 3^t>o, ^) and a standard 
Wiener process {W(t))t>o, we consider the stochastic process in two di- 
mensional space {3^)t>o- 

The system of Ito equations associated to system ffTUl) is given, in the 
equilibrium point P{xo,yo), by 

x{t) = Xo + J^[x{s){hi{x{s)) - h2{x{s))y{s)]ds + gi{x{s),y{s))dW{s), 
y{t) = yo + Jo [{hixis)) - h4{x{s)))y{s) + h{x{s))]ds + g2{x{s), 
y{s))dW{s), 

(12) 

where the first integral is a Riemann integral, the second one is an Ito 
integral and {W{t)}t>o is a Wiener process [17]. 
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The functions gi{x, y) and g2{x, y) are given in the case when we are 
working in the equihbrium state Pg, and they are given by 

gi{x, y) = biix + 612^ + cu, 

(13) 

g2{x, y) = b2lX + 6222/ + C2e, 

where 

Cie = -biiXe - bi2ye, ^ = 1,2, (14) 
and bij G M, i,j = 1, 2. 

3.1 Analysis of Bell model. Numerical simulations. 

3.1.1 Lyapunov exponent method 

Following the algorithm for determining the Lyapunov exponent (An- 
nexe 1) and the description of second order Euler scheme (Annexe 2) in 
Maple 12 software, we get the following results, illustrated in the figures 
below. For numerical simulations we use the following values of param- 
eters: 

ai = 2.5, as = 1, &i = 1, &2 = 0.4, 63 = 0.95, 64 = 2. 
The matrices A and B are given, in the equilibrium point Pi^O, by 

\-b2 + biyi bixi-bsj ' \p a J ' 

with a = a G M, /3 = —2. In a similar way the matrices A and B are 
defined in the other equilibrium point Af^^^-^^^, — )• 




Figure 15: in Pi for ODE ^ Figure 16: {n,x{n,uj)) in Pi for SDE ^ 
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Figure 17: in Pi for ODE ^ Figure 18: {n,y{n,uj)) in Pi for SDK ^ 
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Figure 19: (x(n),y(n)) in Pi for ODE ^ Figure 20: w), w)) in Pi for SDE (O 




Figure 21: in P2 for ODE ^ Figure 22: (n,a;(n,w)) in P2 for SDE (O 
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Figure 23: {n,y{n)) in P2 for ODE ^ Figure 24: {n,y{n,Lu)) in P2 for SDK ^ 

The variation of Lyapunov exponent with the variable parameter 
611 = a is given in Figure 23 for Pi and in Figure 24 for P2. 




Figure 25: (a, A(a)) in Pi for ODE ([T0|) Figure 26: (a, A(a)) in F2 for SDE ^ 



From the figures above, the equihbrium points Pi and P2 are asymp- 
totically stable for all a such that the Lyapunov exponents A (a) < 0, and 
unstable otherwise. So, Pi is asymptotically stable for a G {—00, — 2.02)U 
(1.78, 00) and P2 is asymptotically stable for a G (— 00, — 1.62)U(1.88, 00). 

3.1.2 Lyapunov function method 

For the system of differential equations that describes Bell model, the 
next assertions are true. 

Proposition 3.1 (a) The matrix of the system of differential equations 
that describes the linearized in P2 is 
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A 

where 



ai2 

0-21 ^22 



_ a2(cn&3 - a.2&4) _ ai6i - 0262 a2(&i&4 + &2&3 

^12 — ; ; , 0,21 — , 022 



aihi — 0262 0-2 dibi — a2&2 

(h) If aibi — a2&2 > and 0163 — 0264 > 0, then the equilibrium point P2 
is asymptotically stable. 

□ 

The stochastic model is given using a perturbation around the equi- 
hbrium point P2 (3^2,1/2), in the following way 

dx{t) = x(t){ai — a2y{t))dt + o"i(x(t) — xp^)dW^, 
dy{t) = [(6ix(t) - b-Mt) - b2x{t) + b^]dt + a2{y{t) - yp,)dW\ 

(15) 

with (Ti > 0, (72 > 0. 

The linearized of system ffT^ in (0, 0) is given by 

du{t) = h{u{t))dt + l{u{t))dW{t), (16) 

where u{t) = {ui{t) , U2{t)f , W{t) = {W\t),W'^{t)f and 

h(,,(i\\-( {ai - a2y2)ui{t) - X2a2U2{t) \ 

\{biy2-b2)u,{t) + {b,X2-b,)u2{t)) ' ^''> 

We consider the set = {(t > 0) x M^} and V : D ^ ^ a. function 
of class with respect to t, and of class with respect to the other 
variables. We study the p— exponential stability of the solution (0, 0) of 
the linearized stochastic system ffTUl) . Using Theorem I4.4|, from Anexe 
Al, for the function V : D 

V{t,u) = ]^{ujiu\ + UJ2UI), Ui,U2eM.+ , (19) 
we get the following result. 

Proposition 3.2 // the following relations take place 



qi = Ui{a2y2 - ai- of) > 0, q2 = CU2(&3 - &1Z/2 - crl) > 0, 
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6l7/2 - 62 > 0, t^i = UJ2, 

then 

dV{t,u) = -u{tfQu{t), 

'qi 0' 
,0 q,^ 

The equilibrium point of ( |T5i) asymptotically stable in quadratic 
square (p = 2). 



with Q given by Q 



Proof: From ([HD, ([HD and ([19]), we get 

dV(t u) = ( ("1 ~" "21/2)miW -a;2a2M2 \^ fuJiuA ^ fafuiuj 

= -qiul - q2ul + {uj2{biy2 - ^2) - u^ia;2a2)Miti2- 

If the relations from the proposition take place, then we get 

dV{t,u) = -u{tfQu{t). 

The matrix Q is symmetric and positive defined and has positive 
eigenvalues ri = qi and r2 = q2- Let be = 'nim{gi,g2}- Results 
that 

LVit,u) < -g„||M(t)f. 
and so the equilibrium point is asymptotically stable in square mean. 

□ 

Let us choose the same parameters values for ai, 02, 61, 62? ^3 as on 
the simulation of Lyapunov exponents. We use Maple 12 software for 
the implementation of the second order Euler method. We observe from 
the following graphics that the solution trajectories represent the stable 
characteristic, which validate our theoretical discussion for the system of 
differential equation (fTSll . for the equilibrium point P2. 



13 




0.27 

0.26- 

0.25 

0.24 
0.23 



100 200 300 400 500 600 700 800 



Figure 27: in Pi for ODE ^ Figure 28: (n,a;(n,w)) in Pi for SDE U 




100 200 300 400 500 600 700 800 



Figure 29: in Pi for ODE ^ Figure 30: (n, a;(n, cj)) in Pi for SDE mi) 



4 Conclusions 



In this paper we focused on two important tumor-immune systems, 
presented from stochastic point of view: a Kuznetsov- Taylor model and 
Bell model, that belongs to a general family of tumor-immune stochastic 
systems. We have determined the equilibrium points and we have calcu- 
lated the Lyapunov exponents. A computable algorithm is presented in 
Annexe Al. These exponents help us to decide whether the stochastic 
model is stable or not. We have also proved stochastic stability by con- 
structing a proper Lyapunov function, under a well chosen conditions. 
All our results were also proved using graphical implementation. For 
numerical simulations we have used the second order Euler scheme pre- 
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sented in detail in Annexe A2 and the implementation of this algorithm 
was done in Maple 12. 

In a similar way can be studied other models that derive from model 
given by ( |TOl) . The model given by the SDE ( |T2i) allows the control of 
the system given by the ODE ([T]), with a stochastic process. 

Annexe 

Al Lyapunov exponents and stability in stochastic 2- 
dimensional structures. 

Lyapunov exponent method 

The behavior of a deterministic dynamical system which is disturbed 
by noise may be modeled by a stochastic differential equation (SDE). In 
many practical situations, perturbations are generated by wind, rough 
surfaces or turbulent layers are expressed in terms of white noise, mod- 
eled by Brownian motion. The stochastic stability has been introduced 
in [13] and is characterized by the negativeness of Lyapunov exponents. 
But it is not possible to determine this exponents explicitly. Many nu- 
merical approaches have been proposed, which generally used simulations 
of stochastic trajectories. 

Let {Q, 5", T) a probability space. It is assumed that the cr— algebra 
(5't)t>o such that 

9^, C C 3", Vs < t, s,t e /, 

where I = [0,T], T e (0,oo). 

Let {x{t,u) = {xi{t), X2it j)}t>o be a stochastic process, solution of 
the system of Ito differential equations, formally written as 

dxi{t, u) = fi{x{t, uj))dt + gi{x(t, uj))dW{t, cj), i = 1, 2, (20) 

with initial condition x(0) = xq is interpreted in the sense that 

Xiit.uj) = Xioit.oj) + / fi{x{s,uj))ds + / gi{x{s,uj))dW{s,u), i = 1,2, 
Jo Jo 

(21) 

for almost all u & Q and for each t > 0, where fi{x) is a drift function, 
gi{x) is a diffusion function, fi(x{s))ds, i = 1, 2 is a Riemann integral 
and gi{x{s))dW{s), z = 1,2 is an Ito integral. It is assumed that fi 
and gi, i = 1,2 satisfy thee conditions of existence of solutions for this 
SDE with initial conditions x(0) = ao G M". 
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Let Xq = (xoi, X02) G be a solution of the system 



/i(xo) = 0,2 = 1,2. 
The functions Qi are chosen such that 

gi{xo) = 0, i = 1,2. 
In the following, we will consider 



(22) 



2 



1,2, 



(23) 



where bij G 



The linearized of the system fl^Tl) in xq is given by 



du{t) = Au{t)dt + Bu{t)dW{t), 



where 



u{t) 



Ui{t) 
U2{t) 



A 



Oil ai2 

0'2l 0,22 



B 



bn bi2 
&21 ^22 



dxi 



b - = — 

dxi 



(24) 



(25) 



(26) 



x=xo 



x=xo '^■^J 

The Oseledec multiplicative ergodic theorem [16] asserts the existence 
of two non-random Lyapunov exponents A2 < Ai = A. The top Lyapunov 
exponent is given by 



A = lim sup - log a/ ui{ty + U2{ty- 



By applying the change of coordinates 

ui{t) = r{t) cos^(t), U2it) = r{t) sin^(t), 

for ([21]) and by using the Ito formula for 

1 



(27) 



hi{ui,U2) 
h2{Ui,U2) 



arctan ( — ) =6*, 

Ml 



result the stochastic equations written in the integral form. 



log 



r{t) 
r(0) 



qMs)) + -{qMs)f-q2{e{s)f)ds+ j q2{e{s))dW{s) 

(28 
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e{t) = e{0)+ I {q,{e{s))-q2{e{s))qMs)))ds+ I qMs))dW{s), (29) 
where 



(30) 



qi{0) = ail cos^ 9 + (a^ + a2i) 008 6*8111^ + 022 sin^^, 
92(6') = bii co8^ 9 + {bi2 + 621) C08 9 sin 9 + 622 sin^ 9, 
q3{9) = a2i co8^ 9 + (022 — an) cos 9 sin 9 — ai2 sin^ 9, 
q4:{9) = 621 co8^ 9 + (622 — ^11) C08 9 sin 9 — bi2 sin^ 9. 

As the expectation of the Ito 8tocha8tic integral i8 null, 

E(^l\2{0{s))dW{s))=O, 
the Lyapunov exponent i8 given by 

Applying the 08eledec theorem, if r{t) is ergodic, re8ult8 that 

A = l\qi{9) + ^{q,{9y - q2{9)))p{9)d9, (31) 

where p{9) is the probability den8ity of the proce88 9. 

The probability den8ity i8 the 8olution p{t, 9) of Fokker-Planck equa- 
tion a880ciated to equation 



^^^^'^^ + 4fe(^) - q2{9)q,{9)p{t,9)) - l^(q,(9fp{t,0)) = 0. (32) 



If pit, 9) = p{9), then the stationary solution of (IH^ satisfies the first 
order differential equation 

i-qsiO) + qi{9)qm + q2{9)q,{9))p{9) + ^q^'piO) = Po, (33) 

where p{9) = ^ and 

qm = -{bi2 + 621) sin(2^) - (622 - &11) cos(2^). (34) 
Proposition 4.1 If q4{9) 7^ 0, the solution of equation ( l33l) is given by 
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where K is determined by the normality condition 



2tt 

p{e)d9 = 1, (36) 



and 

D{27t) - 1 
J^^D{u)du 
The function D is given by 



V = ^: . , • (37) 



D{e) = exp{-2 I '?3H-g2Hg4(n) -g4Hg5H ^^)_ (3g) 



□ 

A numerical solution of the phase distribution could be performed by 
a simple backward difference scheme. The function p{6) can be deter- 
mined numerically by using the following algorithm. 

Let us consider N G M+, h = and 

qiii) = ail cos^{ih) + (ai2 + 021) cos{ih) sm{ih) + 022 sm^{ih), 
q2{i) = bii cos^(i/i) + (612 + &21) cos{ih) sm{ih) + 622 sin^(i/i), 
q3{i) = a2i cos^{ih) + (022 — On) cos{ih) sin{ih) — 022 sin^^ih), 
qiii) = &21 cos^(i/i) + (622 — bii) cos{ih) sm{ih) — bu sin^(i/i), 
95(0 = -{bu + &21) sm{2ih) - (622 - &11) cos{2ih), i = 0, 1, iV. 

(39) 

The sequence (p(i))i=o,...,Ar is given by 

p{i) = (p(0) + — )F{t), 

where 



2h{-q3{i) + 92(0^4(0 + 94(095(0) + 94(0^ 
The Lyapunov exponent is A = X{N), where 

^ 1 
A(iV) = Y.^ql(^) + 2(94(0' - 92(0'))P(0/^- 

i=l 

From ( l30l) and ( l35i) we get the following proposition. 
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Proposition 4.2 // the coefficients of the the matrix B are given by 

hi = tt, bi2 = -P, &21 = P, &22 = a, 
then the Lyapunov exponent is given by 

1 11 

A = -(an + 022 + - a^) + -(an - a22)-D2 + -(021 + ai2)-E2, 

where 

D2 = J cos{2e)p{e)de, ^2 = J sm{2e)p{e)de. 



p{e) = Kg{e), K 





1 



tame' 



3(e) = ■^exp{^{{a2i-ai2-aP)9+^{aii-a22) cos(26') + ^(a2i-ai2) sin(26'))). 

□ 



Lyapunov function method 

Let us consider the stochastic system of differential equations given 

by 

dxi{t) = fi{x{t))dt + gi{x{t))dWi{t), % = \X (40) 

where W\, W2 are Wiener processes. Let D = (0, 00) x R^, and V : 
D —>■ a. continuous function with respect to the first component and 
of the class with respect to the other components. Let consider the 
differential operator given by 



2 r^T^/, X ,22 



,,,, , dVit,x) ^^^^dV{t,x) Iv^v^ , ,d'V{t,x) 

1=1 i=l j = l ■' 

(41) 

We suppose that Xq = is an equilibrium point for (HUj) . that means 

/i(0) = ^?.(0) = 0, 1 = 1,2. (42) 

The theorem that gives the conditions for stability of the trivial so- 
lution Xo = in terms of Lyapunov function is given in [T7] . 
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Theorem 4.3 // there is a function V : U ^ and two continuous 
functions u,v : M_|_ IR+ and k > 0, such that for each \\x\\ < k, we 
have 

ui\\x\\)<Vix,t)<vi\\x\\), (43) 

then 

(i) If LV(t,x) < 0, X E (0, /c), then the solution Xq = of (HUjl is stable 

in probability, 

(ii) // there is a continuous function c : IR+ IR+ such that 

LV{t,x)<-ci\\x\\), 
then the solution xq = of fHUl) is asymptotically stable. □ 

In general, the functions /«, Qi, i = 1,2, are nonhnear and the above 
theorem is hard to use. That is why we use the hnearization method for 
the system ( HOj) around the equihbrium point. 

The hnearized system of stochastic differential equation of ( HOl) is 
given by 



r dui{t) = (aiiMi(t) + ai2U2{t))dt + + hi2U2{t))dWi, 

\ du2{t) = (a2iMi(t) + a22U2{t))dt + {b2iUi{t) + 622^2 (^))c^W^2- 

We consider D = {(t > 0) xM^} and V : D a continuous function 
with respect to t and of the class with respect to the other components. 
The theorem that gives the condition that the trivial solution of (jH]) is 
exponential p— stable is given in fl] . 

Theorem 4.4 // the function V : D satisfies the inequalities 

ki\\u\\P < V{t,u) < k2\\u\\P, 

LV{t,u) < -ksWuW, ki>0,p>0, 
then the trivial solution of (1441) is exponentially p— stable for t > 0. □ 

In concrete problems, the next theorem is used. 

Theorem 4.5 // the function V : D satisfies 

(i) LV{u) < 0, then the trivial solution is stable in probability; 

(ii) LV{u) < — c(||u||), where c : M_|_ ]R_|_ is a continuous function, 

then the trivial solution is asymptotically stable; 
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(iii) LV{u) < —q^Qq, where Q is a symmetric matrix, positive defined, 
then the trivial solution is stable in mean square value. 



□ 



For (jH]), the expression of the differential operator LV is given by 



LV{t,u) 



,dV(t,u) , ,dV(t,u 
[aiiUi + ai2U2) — h (a.2iMi + «22M2 



+ 



dui 



+ (&21M1 + 622M2) 



dU2 

2d'Vit,u) 



dul 



A2 The Euler scheme. 

In general 2-dimensional case, the Euler scheme has the form: 

Xi{n + 1) = Xi{n) + fi{x{n))h + gi{x{n))Gi{n), i = 1,2, (45) 
with Wiener process increment 

Gi{n) = Wi{{n + l)h) - Wi{nh), n = 0, N - 1, i = 1,2, 

and Xi{n) = Xi{nh, omega). Gi{n) are generated using Box-Muller method. 

It is shown that the second Euler scheme has the order for weak 
convergence 1, for sufficiently regular drift and diffusion coefficients. 

We assume that fi in fl45p are sufficiently smooth such that the fol- 
lowing schemes are well defined. 

The second order Euler scheme is defined by the relations 



Xi{n + 1) 



d GAnY - h 
Xi{n) + fi{x{n))h + gi{x{n))Gi{n) + gi{x{n)) gi{x{n))— — h 



+ 



dh{x{n)) ^ 1 
dxi{n) 2' 
,dgi{x{n)) ^ 1 
2 



dxi{n) 



+ h{x{n)r-^^^^ + ^g,{x{n)f 



dxi{n) 



dxi{n)dxi{n)\ 2 
d'^gi{x{n)) 1 HGiiri 
dxi{n)dxi{n) 



gi{x{n)) 



dft{x{n)) 
dxiiji) 



1,2, 



where we used the random variables Gi{n), = 1,2. In [12], it is shown 
that these schemes converge weakly with order 2. 
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